{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "影像拼接：\n",
    "1.创建空矩阵（相当于底图）\n",
    "2.将待镶嵌的影像读入，并添加到空矩阵中"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_img_1=\n",
    "path_img_2=\n",
    "dset_subs1=gdal.Open(path_img_1)\n",
    "dset_subs2=gdal.Open(path_img_2)\n",
    "geo_trans_subs1=dset_subs1.GetGeoTransform()\n",
    "geo_trans_subs2=dset_subs2.GetGeoTransform()\n",
    "print(geo_trans_subs1)\n",
    "print(geo_trans_subs2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### image 1\n",
    "x_max_subs1=geo_trans_subs1[0]+dset_subs1.RasterXSize*geo_trans_subs1[1]\n",
    "y_min_subs1=geo_trans_subs1[3]+dset_subs1.RasterYSize*geo_trans_subs1[5]\n",
    "print(x_max_subs1,y_min_subs1)\n",
    "extent_subs1=[geo_trans_subs1[0],x_max_subs1,y_min_subs1,geo_trans_subs1[3]]\n",
    "print(extent_subs1)\n",
    "\n",
    "### image 2\n",
    "x_max_subs2=geo_trans_subs2[0]+dset_subs2.RasterXSize*geo_trans_subs2[1]\n",
    "y_min_subs2=geo_trans_subs2[3]+dset_subs2.RasterYSize*geo_trans_subs2[5]\n",
    "print(x_max_subs2,y_min_subs2)\n",
    "extent_subs2=[geo_trans_subs2[0],x_max_subs2,y_min_subs2,geo_trans_subs2[3]]\n",
    "print(extent_subs2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min_mosaic=min(geo_trans_subs1[0],geo_trans_subs2[0])\n",
    "x_max_mosaic=max(x_max_subs1,x_max_subs2)\n",
    "y_min_mosaic=min(y_min_subs1,y_min_subs2)\n",
    "y_max_mosaic=max(geo_trans_subs1[3],geo_trans_subs2[3])\n",
    "extent_mosaic=[x_min_mosaic,x_max_mosaic,y_min_mosaic,y_max_mosaic]\n",
    "print(extent_mosaic)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对空间分辨率进行统一"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_res_mosaic,y_res_mosaic=geo_trans_subs1[1],geo_trans_subs1[5]\n",
    "print(x_res_mosaic,y_res_mosaic)\n",
    "path_img_2_resam=\n",
    "warp_resam=gdal.Warp(destNameOrDestDS=path_img_2_resam,srcDSOrSrcDSTab=path_img_2,xRes=20,yRes=20,)\n",
    "dset_subs2_resam=gdal.Open(path_img_2_resam)\n",
    "geo_trans_subs2_=dset_subs2_resam.GetGeoTransform()\n",
    "print(geo_trans_subs2_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_size_mosaic=(x_max_mosaic-x_min_mosaic)/x_res_mosaic\n",
    "y_size_mosaic=(y_max_mosaic-y_min_mosaic)/y_res_mosaic\n",
    "print(x_size_mosaic,y_size_mosaic)\n",
    "x_size_mosaic,y_size_mosaic=round(x_size_mosaic),round(y_size_mosaic)#取整\n",
    "x_max_mosaic=x_min_mosaic+x_size_mosaic*x_size_mosaic\n",
    "y_min_mosaic=y_max_mosaic+y_size_mosaic*y_size_mosaic\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geotrans_mosaic=[x_min_mosaic,x_res_mosaic,0,y_max_mosaic,0,y_res_mosaic]\n",
    "print(geotrans_mosaic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#定义矩阵，三个维度分别是波段数、行数和列数\n",
    "img_arr_mosaic=np.zeros(shape=(dset_subs1.RasterCount,y_size_mosaic,x_size_mosaic))\n",
    "print(img_arr_mosaic)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_start_subs=int((geo_trans_subs1[3]-y_max_mosaic)/y_res_mosaic)\n",
    "col_start_subs=int((geo_trans_subs1[0]-x_min_mosaic)/x_res_mosaic)\n",
    "print(row_start_subs,col_start_subs)\n",
    "row_start_subs2=int((geo_trans_subs2[3]-y_max_mosaic)/y_res_mosaic)\n",
    "col_start_subs2=int((geo_trans_subs2[0]-x_min_mosaic)/x_res_mosaic)\n",
    "print(row_start_subs2,col_start_subs2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#将影像1的矩阵赋值到空数组中\n",
    "img_arr_mosaic[:,row_start_subs:row_start_subs+dset_subs1.RasterYSize,\n",
    "               col_start_subs:col_start_subs+dset_subs1.RasterXSize\n",
    "               ]=dset_subs1.ReadAsArray()\n",
    "#将影像2的矩阵赋值到空数组中\n",
    "img_arr_mosaic[:,row_start_subs2:row_start_subs2+dset_subs2.RasterYSize,\n",
    "               col_start_subs2:col_start_subs2+dset_subs2.RasterXSize\n",
    "               ]=dset_subs2.ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_mosaic=\n",
    "driver=gdal.GetDriverByName('GTiff')\n",
    "dset_mosaic=driver.Create(path_mosaic,\n",
    "                          xsize=x_size_mosaic,\n",
    "                          ysize=y_size_mosaic,\n",
    "                          bands=dset_subs1.RasterCount,\n",
    "                          eType=gdal.GDT_Int16)\n",
    "dset_mosaic.SetGeoTransform(geotrans_mosaic)\n",
    "dset_mosaic.SetProjection(dset_subs1.GetProjection())\n",
    "\n",
    "for i in range(dset_subs1.RasterCount):\n",
    "    outband=dset_mosaic.GetRasterBand(i+1)\n",
    "    outband.WriteArray(img_arr_mosaic[i])\n",
    "dset_mosaic=None"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
